Linear density response function in the projector-augmented wave method: 
Apphcations to sohds, surfaces, and interfaces 
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We present an implementation of the linear density response function within the projector- 
augmented wave (PAW) method with applications to the linear optical and dielectric properties 
of both solids, surfaces, and interfaces. The response function is represented in plane waves while 
the single-particle eigenstates can be expanded on a real space grid or in atomic orbital basis for 
increased efficiency. The exchange-correlation kernel is treated at the level of the adiabatic local 
density approximation (ALDA) and crystal local field effects are included. The calculated static 
and dynamical dielectric functions of Si, C, SiC, AlP and GaAs compare well with previous cal- 
culations. While optical properties of semiconductors, in particular excitonic effects, are generally 
not well described by ALDA, we obtain excellent agreement with experiments for the surface loss 
function of the Mg(OOOl) surface with plasmon energies deviating by less than 0.2 eV. Finally, we 
apply the method to study the influence of substrates on the plasmon excitations in graphene. On 
SiC(OOOl), the long wavelength tt plasmons are significantly damped although their energies remain 
almost unaltered. On Al(lll) the vr plasmon is completely quenched due to the coupling to the 
metal surface plasmon. 

PACS numbers: 73.20.Mf, 71.15.-m, 78.20.-e. 
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I. INTRODUCTION 

Time-dependent density functional theory (TDDFT)^ 
has been widely used to calculate optical excitations in 
molecules and clusters as well as the optical and electron 
energy loss spectra of bulk semiconductors, metals and 
their surfaces^. The excitation energies and oscillator 
strengths of both single-particle and collective electronic 
excitations are determined by the frequency-dependent 
linear density response function x(r, r', w) giving the den- 
sity response at point r to first order in a time-dependent 
perturbation of frequency w applied at point r', 



5n{Y,uj) = / (irx(r,r',i:j)(5V'oxt(r',cj). 



(1) 



For finite systems, x can be efficiently calculated by in- 
verting an effective Hamiltonian in the space of particle- 
hole transitions. For the practically relevant case of 
frequency-independent exchange-correlation kernels this 
formulation leads to the well known Casida equation-. 
For extended systems, it is more convinient to express x 
in a basis of plane waves^^— where it has the generic form 
XGG'(qiw), with G being reciprocal lattice vectors and 
q being wavevectors in the first Brillouin zone (BZ). 

In this paper we focus on the electronic response func- 
tion of extended systems treating electron-electron in- 
teractions at the level of the random phase approxi- 
mation (RPA) and the adiabatic local density approx- 
imation (ALDA). For many extended systems such a 
description is insufficient to account for optical excita- 
tions because the electron-hole attraction is not prop- 
erly accounted for. However, dielectric properties, in 
particular collective plasmon excitations, are generally 
accurately reproduced by this approach^*- , and quanti- 



tative agreement with electron energy loss experiments 
have been reported for bulk metals^ii^, surfaces^ii^^, 
graphene-based systemsi^iii^, semiconductorsi^ii^ and 
even supercondutorsil. Furthermore, the accurate eval- 
uation of the density response function at the RPA or 
ALDA level is a prerequisite for implementation of most 
post-DFT schemes, such as RPA correlation energy—, 
exact-exchange optimized-effective-potential methoda^^, 
the GW approximation for quasi-particle excitationsS2i2i, 
and the Bethe-Salpeter equation^-:^^ for optical excita- 
tions. 

Here we present an implementation of the density 
response function within the electronic structure code 
GPAW^'^'^'' which is based on the projector augmented 
wave (PAW) methodology^^'^® and represents wave func- 
tions on real space grids or in terms of linear combina- 
tions of atomic orbitals (LCAO)^'^. Within the PAW for- 
malism one works implicitly with the all-electron wave 
functions and has access to the (frozen) core states. This 
makes the method applicable to a very broad range of sys- 
tems including materials with strongly localized d or f 
electrons which can be problematic to describe with pseu- 
dopotentials. An additional advantage of the PAW for- 
malism, with respect to linear response theory, is that the 
optical transition operator in the long wavelength limit 
can be obtained directly due to the use of all-electron 
wavefunctions^^. The non-interacting response function, 
X°, is built from the single-particle eigenstates obtained 
either on a real space grid, which is the standard repre- 
sentation in the GPAW code, or in terms of a localized 
atomic orbital (LCAO) basis. We have found that the 
latter choice reduces the computational cost of con- 
siderably while still preserving the high accuracy of the 
grid calculation. 
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The method is used to calculate the macroscopic di- 
electric constants of a number of bulk semiconductors, 
showing very good agreement with previous calculations 
as well as experiments. For the surface plasmons of the 
Mg(OOOl) surface we find, in agreement with previous 
studies, that the ALDA kernel lowers the plasmon en- 
ergies by around 0.3 eV relative to the RPA values and 
thereby reduces the deviation from experiments from 4% 
to 1-2%. Very good agreement with experiments is also 
found for the plasmon energies of graphene which are 
shown to exhibit a linear dispersion with a value of 4.9 eV 
in the long wave length limit. The deposition of graphene 
on a SiC substrate is shown to have little effects on the 
plasmon energies but leads to significant broadening of 
the plasmon resonances. In contrast deposition on an Al 
surface completely quenches the graphene plasmons due 
to strong non-local electronic screening. 

The rest of this paper is organized as follows. Section 
II introduces the theoretical framework, where the PAW 
methodology, the density response function for both fi- 
nite q and q — )■ 0, and the ALDA kernel in the PAW 
method are discussed. The details of the implementation 
and parallelization in GPAW and other technical details 
are presented in section III. Section IV presents appli- 
cations for optical properties and plasmon excitations of 
bulk and surfaces, where comparison with other calcula- 
tions and experiments are given. Our recent investigation 
on the effect of a semiconducting and metallic substrate 
on the plasmon excitations in graphene is also briefly 
discussed in this section. Finally, a summary is given in 
section V. 



II. METHOD 

A. Basics of the PAW formalism 

In the PAW formalis m^^i^^ , a true all-electron Kohn- 
Sham wavefunction ipnk is obtained by a linear trans- 
formation from a smooth pesudo-wave-function ipnk via 
ipnk — Ti'nk- The transformation operator is chosen in 
such a way that the all-electron wavefunction ipnk is the 
sum of the pseudo one V'nk and an additive contribution 
centered around each atom written as 

^nk(r) = ^„k(r) + (p^l^nk) - Ra) - C(r - R,)] 

a,i 

(2) 

The pseudo-wave-function ipnk matches the all-electron 
one V'nk outside the augmentation spheres centered on 
each atom a at position R^. Their differences inside the 
augmentation region are expanded on atom-centered all- 
electron partial waves (j)f and the smooth counterparts 
(j)f. The expansion coefficient is given by (p^lV'nk), where 
is chosen as a dual basis to the pseudo-partial wave 
and is called a projector function. A frequently occuring 
term is the all-electron expectation value for a semilocal 



operator A written as 

(V'nkl^lV'nk) = ('0nk|^|V'«k) 

+ 5^(^„kipn(p?i^nk)[(0?i^i0p - mAi^'^)] (3) 

a,ij 

B. Density response function and dielectric matrix 

A key concept in TDDFT is the density response func- 
tion It is defined as x(r,r',w) = 6n{r , uj) / SVcxti'r' , io) , 
where T4xt is the external perturbing potential and 6n is 
the induced density under the perturbation. For periodic 
systems, x can be written in the form 

^ BZ 

^ q GG' 

(4) 

where G,G' are reciprocal lattice vectors, q is a wave 
vector restricted to the first BrouUion Zone (BZ), Nq is 
the number of q vectors and CI is the volume of the real 
space primitive cell. 

The density response function of the interacting elec- 
tron system, Xi can be obtained from the non-interacting 
density response function of the Kohn-Sham system, x", 
and a kernel, K , describing the electron-electron interac- 
tions by solving a Dyson-like equation 

XGG'{(l,Uj) = XGG'i^i,^) 

+ E XGGi(q,'^)^GiG2(q)xG2G'(q,w). (5) 

G1G2 

The expression for the non-interacting density re- 
sponse function in the Bloch representation of Adler and 
Wiser— 1^, is 

2 

XgG' (q, ^) = qYI (^"k - /«'k+q) 
k,7in' 

?lnk,n'k+q(G)n*^ (G') 

X ■ , (6) 

W + e„k - En'k-l-q + IT] 

where 

«„k,„'k+q(G) = (V^nkle-'^^+^^-n^n'k+q) (7) 

is defined as the charge density matrix. Its evaluation 
within the PAW formalism is explained in detail in the 
following subsection. e„k, fnk and ipnk are the Kohn- 
Sham eigen-energy, occupation and wave function for 
band index n and wave vector k, and rj is a broadening 
parameter. The summation over k runs all over the BZ 
and fnk = 1 is satisfied for the occupied states. The 
factor of 2 accounts for spin (we assume a spin-degenerate 
system). 

The kernel in Eq. (O consists of both a Coulomb and 
an exchange-correlation(xc) part. The Coulomb kernel is 
diagonal in the Bloch representation and written as 

n / > 47r „ , ^ 

^SxG.(q)= |^q-^^G,G., (8) 
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while the xc kernel evaluated within ALDA is given by 
1 



with 



-j(Gi-G2)-r 



/xc[ri(r)] = 



, (9) 



(10) 



no(r) 



Details on the evaluation of the xc kernel in the PAW 
method can be found in a following subsection. 

The Fourier transform of the microscopic dielectric ma- 
trix, defined as e~^(r,r',a;) = (5Vtot(r,a;)/JT4xt(r', w), is 
related to the density response function via 



47r 



G|2 



XGG'(q,t^) (11) 



where x is obtained from according to Eq. ([5]). The 
off-diagonal elements of the Xgg' nis-trix describes the 
response of the electrons at wave vectors different from 
the external perturbing field and thus contain informa- 
tion about the inhomogeneity of the microscopic response 
of electrons known as the 'local field effect—. The macro- 
scopic dielectric function is defined as 



eji/(q,'^) 



1 



](q,^) 



(12) 



■=00 



and is directly related to many experimental properties. 
For example, the optical absorption spectrum (ABS) is 
given by ImeM(q — >■ 0, w). The electron energy loss spec- 
trum (EELS2S) is propotional to — Im(l/eM)- Both spec- 
tra reveal information about the elementary electronic 
excitations of the system. EELS is especially useful in 
probing the collective electronic excitations, known as 
plasmons, of bulk and low-dimensional systemaSi. 



C. Charge density matrix in the PAW method 

In this subsection, we will discuss the charge density 
matrix n„k,n'k+q(G), which is defined in Eq. ([7]) and is 
a crucial quantity for the evaluation of ■ Care must 
be taken for the long wavelength limit (q -> 0) since the 
Coulomb kernel, 47r/|q-|- Gp, diverges at q — >■ and 
G = 0; while the charge density matrix approaches zero 
at this limit. As a result, we separate the discussion into 
two parts: finite q and q — >■ 0. 



1. Finite q 

Considering the transformation between the pseudo- 
wavefunction and the all-clcctron wavefunction in Eq. ([2]) 
and employing Eq. ([3]) yields 

«nk,ri'k+q(G) = ^ink,ri'k+q(G) (13) 

+ E(^»k|pn(p,1^n'k+q)Q°,(q + G) 



with 

n„k,„'k+q(G) = (7^„k|e-'('i+°)-"|^n'k+q) (14) 
Qf^{K) ^ (C|e-^^-"-|0,") - (^,?|e-K--|0,") (15) 

and K = q + G. 

The pseudo-density matrix in Eq. ((T^ is calculated 
using a mixed space scheme. First, the cell periodic func- 
tion V^[^(i")'0n'k+q(r)e^*''''' is evaluated on a real-space 
grid; then it is Fourier transformed to get 



nk,n'k+q(G) = T Ck(r)i^n'k+q(r)e 



(16) 



The augmentation part in Eq. (jlSp is calculated on 
fine one-dimensional radial grids centered on each atom. 
Such fine grids are required to represent accurately the 
oscillating nature of the all-electron partial wave in the 
augmentation region. The plane wave term e 
panded using real spherical harmonics by 



-iK r 



-iKr 



An 



Im 



-iyjii\K\r)Yi„,{i)YUK), (17) 



where ji is spherical Bessel function for angular momen- 
tum I and K — K/|K|. Combining the above equa- 
tions and the expression for the partial wave = 
<l>n i i'^)^hmi{r), we can write 



dr 



(18) 



2. Long wave length limit 



In the long wave length limit, the G ^ components 
of the density matrix ?T.„k,n'k+q(G) remain the same as 
that for finite q. Only the G = components need to be 
modified and are written as 



n„k,n'k+q(0)|q^O = (l/lnkje *'*''"|'0n'k+q)q^O- 



(19) 



In Ref. [SCj, the above so called longitudinal form is de- 
rived in the PAW framework by using Taylor expansion 
of the e*'! "" to the first order. Here we adopt an alter- 
native but equivalent form which can be derived using 
the second order k ■ p perturbation theorjsSi as described 
below. 

Expressing the wavefunction using Bloch's theorem as 
'0nk(r) = u„k(r)e*'' '", where M„k(r) is the periodic Bloch 
wave, the dipolc transition element in Eq. (|19p becomes 



(V'nkle I V'n'k+q) = (Unk|Un'k+q) 



(20) 



For vanishing q, the wavefunction for |Mn'k-i-q) can be 
obtained in terms of those for \umk) through second order 
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perturbation theory: 



\u 



n'k+q/ 



\U„'k) 



2^ — : l"mk) (21) 



^Ti'k ^mk 



The perturbing potential y in the above equation is ob- 
tained through 



V = H{k + q) - H{k) = -iq • (V + ik) 



where 



ff(k) 



(V + zk)2 + y(r) 



(22) 



(23) 



is the k ■ p hamiltonian2i and V{r) is the effective Kohn- 
Sham potential. 

Combining Eq. (|20p - ([2^ . the charge density matrix 
at the long wavelength limit becomes 



^nk.n'k+q (0)|q^0 



-^q- (V^nk|V|V^n^k) 



(24) 



The above expression for the charge density matrix in the 
PAW method has an advantage over the pseudopotential 
method, where the nabla operator has to be corrected by 
the commutator of the non-local part of pseudopotential 
with the position operator r^*. In the PAW method, the 
matrix element (^Ankl Vji/in/k) is given by 



(^AnklVIV'n'k) = (V'«k|V|l/j„'k) 



(25) 



In GPAW, where the pseudo wave functions, 4'nk, 
represented on a real space grid, the first matrix element 
is calculated using a finite difference approximation for 
the nabla operator. The augmentation part is evaluated 
on fine one dimensional radial grids. The nabla operator 
combined with partial waves (/)f(r) = 4'n^i-^{^)Yiimi{i^) 
and 0°(r) = (l>n^i^{r)Yi^,nA^) is written as 



d 



dr 



dr Qy 



f^2 



(26) 



Since real spherical harmonics are employed, we get 

ar_ xyz /4^ 
CIT r r r \ 6 

Substitute the above equation into Eq. ([2^]) and split the 
integration into radial and angular parts,we get for the 



x-component 

+ Jdr r^C.h^ / dvYi,„,y-^^^y^Yi,„^.) 

(28) 

The derivation for the y- and z-component and for the 
pseudo-partial-wave follows in a similar way. 

D. The ALDA xc kernel in the PAW method 

The ALDA xc kernel, expressed in Eq. © , is evaluated 
using the all-electron density, which takes the form 

n{r) = h{r) + ^K(r - R„) - n'^(r - R„)], (29) 

a 

where 



n(r) 



E 



/„k|Vi„k(r)P+^n^(|r-RJ), (30) 



n'^(r) = l^Dt,<t>tm{r)+n^cir), 
n^{r) = 5]i^^,.0^(r)0^"(r)-t-n^(r), 



(31) 
(32) 



with Df^ = E„k(^"k|pf)/„k(p,"|V^„k). Here n^(r) is the 
all-electron core density and h^{r) can be chosen as any 
smooth continuation of n'^(j) inside the augmentation 
sphere since it will be canceled out in Eq. (PU]). 

The ALDA xc kernel can also be separated into smooth 
and atom-centered contributions 

;.xc— ALDA 



K. 



-ALDA 



G1G2 



^xc— ALDA 
G1G2 



E^^G^G. 



(33) 



The smooth part is constructed from pseudo-density and 
by utilizing a Fourier transform 



K: 



xc-ALDA 
G1G2 



1 



(34) 



{UMr)]}\G,-a, 
The atom-centered contribution is evaluated on ID grids 

1 



AK. 



G1G2 



x[/xcK]-/xc[n"] 



(35) 



III. NUMERICAL DETAILS 



In this section we describe the most important numer- 
ical and technical aspects of our implementation; in par- 
ticular the Hilbert transform used to obtain x° from the 
dynamic form factor (spectral function) and the applied 
parallelization scheme. 
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A. Symmetry 

For each wave vector q, the evaluation of involves a 
summation over occupied and empty states in the entire 
BZ. By exploiting the crystal symmetries, however, we 
need only calculate the wave functions and energies in 
the irreducible BZ. This is because the wave function at 
a general /c-point can always be obtained from a wave 
function in the irreducible part of BZ by application of 
a symmetry transformation, T. In general we have the 
relation 

V'«,Tk(r) = V^„,k(r-ir) (36) 

where k belongs to the IBZ. The above relation can be 
directly verified by considering how the right hand side 
transforms under lattice translations. In addition to the 
crystal symmetries, time reversal symmetry applies to 
any system in the absence of magnetic fields 

^_k(r) = ^i:(r) (37) 

B. Hilbert transform 

Rather than constructing x'^ directly from Eq. ([6|) we 
obtain it as a Hilbert transform of the (non-interacting) 
dynamic form factor, 5'0| 33,34 r^j^^ latter is given by 

2 

k,nn' 

X «nk,n'k+q(G)n*k,„'k+q(Gl')- (38) 

In practice S^{uj) is evaluated on a uniform frequency 
grid extending from to around 40-60 eV with a grid 
spacing in the range 0.01-0.1 eV, and the delta functions 
are approximated by triangular functions following Ref. 
|32| . The non-interacting response function is obtained as 

/•OO 

Jo 



Ll) — Lu' + it] uj + uj' + irj 

The above Hilbert transform is performed directly on the 
frequency grid setting the broadening parameter rj equal 
to the grid spacing. 

C. LCAO vs grid calculations 

It is well known that the use of localized atomic or- 
bitals as basis functions can significantly reduce the com- 
putational effort of groundstate electronic structure cal- 
culations. For calculations of the density response func- 
tion the use of localized basis functions is complicated 
by the fact such basis sets are typically not closed un- 
der multiplication^^"=2I. As a consequence the size of the 
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FIG. 1: (Color online) The imaginary part of the dielec- 
tric function (a) and energy loss function (b) of graphene at 
q = 0.046 along T — M direction of its surface BrouUion 
zone (SBZ) calculated with 3D uniform grid (GRID, black 
solid line) and localized atomic orbital (LCAO) using dzp 
(red dashed line) and qztp (blue dash-dotted line) basis, re- 
spectively . 



product basis needed to represent the response function 
grows as iV^ , where is the number of basis functions 
used to represent the wave functions (we note that for 
strictly localized basis functions, the effective size of the 
"product basis" grows only linearly with the system size 
because pair densities of non-overlapping orbitals van- 
ishes, however, the prefactor is typically very large). A 
further challenge is the computation of the Coulomb in- 
teraction kernel, l/|r — r'|, in the product basis leading 
to six-dimensional multi center integrals. These inter- 
grals must be performed cither by using efficient Poisson 
solvers or by resorting to analytical techniques. The lat- 
ter is extensively used in quantum chemistry codes ap- 
plying Gaussian basis sets. 

For these reasons we have chosen to represent the 
density response function in a plane wave basis. The 
plane wave basis is closed under multiplication and the 
Coulomb kernel is simply given by Eq. ([5]). However, 
we still keep the advantage of using an LCAO as basis 
in the calculation of the Kohn-Sham wave functions and 
energies which enter the construction of x*^ ^ Apart from 
reducing the computational effort of the groundstate cal- 
culation (which must include many unoccupied bands). 
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FIG. 2: (Color online) The energy loss function of Mg(OOOl) 
surface at q — 0.07 A^^ along F — M direction of its SBZ 
calculated with 3D uniform grid (GRID, black solid line) and 
localized atomic orbital (LCAO) using dzp (red dashed line) 
basis. The dzp basis used here includes double-zeta orbitals 
of 3s and 3p atomic orbitals as well as one d-type Gaussian 
polarization function. 



the storage requirements for wave functions become much 
less than for corresponding grid or plane wave calcula- 
tions. This is because the LCAO coefficients provide a 
more compact representation of the wave functions, in 
particular for open structures containing large vacuum 
regions, and because significantly fewer unoccupied wave 
functions result from the LCAO calculation (for a fixed 
energy cut-off). 

Compared to plane waves or real space grids, LCAO 
calculations employing standard basis sets usually give 
a less accurate but often acceptable description of the 
occupied and low-lying unoccupied wave functions and 
energies. For higher-lying unoccupied states, blue shifts 
are expected due to the (unphysical)confinement imposed 
by the localized basis set, and the continuum is broken 
into discrete bands. Despite these effects, we have found 
that the use of LCAO wave functions instead of grid wave 
functions has rather little effect on the dielectric function 
- at least in the relevant low energy regime. 

As an example Fig. [T] shows the absorption spectrum 
(a) and EELS spectrum (b) of graphene calculated us- 
ing wave functions and energies from a grid calculation 
and from an LCAO with double-zeta polarized (dzp) and 
quadruple zeta triple polarized (qztp) basis, respectively. 
The unit cell is the primitive cell of graphene containing 
two carbon atoms and with 20 A vacuum. The BZ is sam- 
pled on a 64 X 64 Monkhorst-Pack grid. The number of 
bands included are 60 for the grid and LCAO (qztp) basis 
and 26 for the LCAO (dzp) basis. In all three cases this 
corresponds to inclusion of states with energy below 40 
eV. The response function is evaluated at the RPA level 
including local field effects up to a plane wave cut-off of 
150 eV. 



For excitation energies below 10 eV, the LCAO results 
agree remarkably well with the grid calculations. The 
TT — > vr* absorption peak at around 4 eV in panel (a) 
and the tt plasmon around 5 eV in panel (b) are well 
reproduced in LCAO calculations. For energies above 
10 eV, we observe slight deviations, however, the overall 
agreement is remarkable for the entire energy range. In 
particular, the cr — > cr* transition at around 14 eV in 
panel (a) and the a plasmon around 17 eV in panel (b) 
are clearly visible, although in the LCAO calculation the 
latter is splitted into two peaks. 

Fig. [5] shows another example of the Mg(OOOl) sur- 
face, which is modeled by a slab of 16 layer. The energy 
loss function calculated with 3D grids is characterized by 
two peaks at around 7.5 and 11 eV, which correspond to 
the surface and bulk plasmons, respectively. Again the 
LCAO (dzp) calculation reproduces the grid results quite 
accurately, the only discrepancy being the slight discrep- 
ancy (0.1 eV) of the peak just below 10 eV. Note that 
dzp basis used in this case includes double-zeta orbitals of 
3s and 3p atomic orbitals as well as one d-type Gaussian 
polarization function. The inclusion of the d-type orbital 
in the basis function is crucial for the correct description 
of both the electronic structure (fx, band structure and 
density of states) and the surface plasmon of the Mg sur- 
face. The response function calculation presented here is 
performed at the RPA level, with summation over bands 
up to 15 eV and inclusion of local fields up to 50 eV 
plane wave cutoff. The frequency grid spacing is 0.1 eV. 
The response function is not fully converged with these 
parameters. Please refer to the next section for the con- 
verged results for Mg(OOOl) surface. 

D. Storage of wave functions 

For ground state calculations performed using grid 
based wave functions, the entire set of occupied and un- 
occupied wave functions might be too large to be stored 
on disk, making the separation of the ground state and 
response function calculations impossible. In this case, 
the response function, or more precisely, the dynamical 
form factor of Eq. (pS)) . is constructed as the wave func- 
tions are calculated. 

In the LCAO mode, only the expansion coefficients of 
the wave functions in terms of the localized basis func- 
tions are calculated and stored. Since this representation 
is significantly more compact than the grid representa- 
tion, the entire set of wave functions can be calculated 
and stored at once, and the calculation of the response 
function can be performed as a post-processing step. 

E. Parallelization 

The calculation of the response function involves 
the three steps: evaluation of the spectral function 
5QQ/(q, oj) according to Eq. (|38p . Hilbert transform fol- 
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FIG. 3: Schematic illustration of the applied parallelization 
scheme. Each box represents a single CPU. (a) The calcula- 
tion of S'QQ/(q, a;) is performed in parallel over wave vectors, 
k, (or bands, n, for large cells) and frequencies u. (b) The 
Hilbert transform is parallelized over G. (c) Finally the Dyson 
equation is solved by parallelizing over the frequencies. 



lowing Eq. ([M)) and solving Dyson's equation Eq. ([SJ. 
Fig. [3] illustrates the parallelization scheme applied for 
each of these three steps. 

It is natural to parallelize the evaluation of "S'qq/ (q, w') 
over k-points (or bands for few k-point calculations). On 
the other hand, the size of the matrix is often too large 
to be handled on a single CPU. In such cases each CPU 
only calculates on a part of the frequency grid. This 
leads to the two-dimensional parallelization scheme il- 
lustrated in Fig. [Sja). Finally the full S'qq, (q, w') is 
obtained by summing over fc-points, i.e. summing up 
the columns in Fig. |3^a). Since the Hilbert transform 
involves a frequency convolution it is convinient to redis- 
tribute the data from parallelization over uj to over G. 
Finally, the Dyson equation is done separately for each 
frequency point and is therefore parallelized over w, as 
shown in panel (c). 



IV. RESULTS 



TABLE I: The static macroscopic dielectric constants e cal- 
culated using the PAW method on the RPA level with- 
out local field (NLF) and including local field (LF) efi^ect. 
These values are compared with other PAW calculations^ 
and experimental^. 
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15.12 


13.67 
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FIG. 4: Imaginary part of the dynamical dielectric function 
of bulk silicon. The arrows indicate the absorption onset and 
the position of main and secondary peaks, respectively, as 
extracted from Ref. l30l. 



In this section, the density response function method is 
applied to study the optical properties and plasmon ex- 
citations of solids. They are usually measured by optical 
and electron energy loss spectroscopy (EELS), which are 
related to Imej\/ and — Im[l/eA/], respectively. For ex- 
tended systems, the two kinds of spectroscopy give quite 
distinct spectra. The optical absorption spectrum (ABS) 
is determined by single-particle excitations while EELS is 
dominated by collective electronic excitations, plasmons, 
which are defined as em — ^ 0. 



A. Optical properties 

Table I shows the calculated RPA static dielectric func- 
tion in the optical limit for five semiconductors (C, Si, 
SiC, AlP, GaAs). We use the same lattice constants as 
in Ref. |30 and a grid spacing of 0.2 A. A Monkhorst-Pack 
grid of 12 X 12 X 12 and 60 unoccupied bands are used. 
We use a Fermi temperature of 0.001 eV in the ground 
state LDA calculation and a broadening parameter {rj) 



of 0.0001 eV in x*^. Note that in this case we calculate 
the static response function directly from Eq. ([5]), i.e. 
we do not use the Hilbert transform. For calculations 
including local field effects, a cutoff of 150 eV is used. 
The dielectric constants obtained both with and with- 
out local fields agree to within 0.1 with previous PAW 
calculations^. The only exception is GaAs for which 
our dielectric constant is 0.4 larger. This deviation could 
come from differences in the PAW setups for Ga or As. 
The inclusion of local fields lowers the dielectric constant 
by 10-15% in agreement with earlier reports^S. These ob- 
tained values are, however, generally larger than the ex- 
perimental values due to the underestimated band gaps 
by LDA. Inclusion of the ALDA kernel only increases the 
dielectric constant further, and are not reported here. 

Fig. m shows the dynamical dielectric function for 
Si. Compared to the calculations for the static dielec- 
tric constant, a significantly denser k-point sampling of 
80 X 80 X 80 is employed here to resolve the finer details 
in the spectrum. A total of 36 unoccupied bands are 
use in the construction of x°- Local field effects are not 
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included and rj is set to 0.01 eV. The onset of absorp- 
tion and the position of the two characteristic peaks in 
the absorption spectrum compare very well with previous 
RPA calculations'^^ as shown by the arrows in the figure. 
However, it is quite different from the experimental ab- 
sorption spectrum'^^ which exhibits an absorption onset 
at ~ 0.5 eV larger than predicted by our calculation, and 
shows a double peak around 3.8 eV. The disagreement 
with the experimental spectrum is due to the underesti- 
mation of the band gap by LDA and the fact that RPA 
does not include the electron-hole interaction. 



B. Plasmon excitations 

In contrast to the optical excitations, like the Si ab- 
sorption spectrum discussed in the previous section, plas- 
mon excitations are generally well described by RPA and 
TDLDA. In the first part of this subsection, a classi- 
cal calculation of the surface plasmons of a Mg(OOOl) 
surface--^ is reproduced and the results are in good agree- 
ment with previous reports^^. In the second part, our re- 
cent investigation of the effect of substrates on the plas- 
monic excitations of graphene is summarized. 

Plasmon excitations appear as strong peaks in the elec- 
tron energy loss spectrum (EELS) which is directly re- 
lated to the imaginary inverse dielectric function, 

-Ime (q,w) = ---^ImxG=o,G'=o q,a;) (40) 

|q| 

For excitations at surfaces, a surface loss function can be 
defined as^^ 

g{q,uj) = J I dzdz'xG||=G;|=o(2;,z';q,w)elil(^+^') 

where || and z correspond to directions parallel 
and perpendicular to the surface, respectively, and 
XG|| g|| {z, z'\ q, uj) is the Fourier transform of xgg' (q, ^) 
in the z-direction. 
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FIG. 5: (Color online) Surface loss function of the Mg(OOOl) 
surface along the T — M direction of the surface BZ calculated 
using RPA (a) and TDLDA (b). In both cases |q| increases 
from bottom to top. (c) Surface plasmon dispersion for both 
the T — M and T — K directions. Results from this work (filled 
dots) compare well with other calculations (hollow dots^) and 
experiments^. 



C. Surface plasmons of Mg(OOOl) 

Fig. [S]shows the surface loss function of the Mg(OOOl) 
surface along the T — M direction of the surface BZ calcu- 
lated within RPA (panel a) and TDLDA (panel b). The 
Mg surface is modeled by a slab of 16 layers as in previous 
calculations'*- , and a vacuum region of 40 A. Such thick 
slab and vacuum region is necessary to avoid splitting of 
the surface plasmon peak due to coupling between the 
surface plasmons at the two sides of the slab. The LDA 
wave functions are calculated on a uniform grid with a 
grid spacing of . 24 A and a 64 x 64 x 1 Monkhorst-Pack k- 
point sampling. For the response function calculations we 
include 200 bands (including 16 occupied bands) and use 
a broadening parameter of 0.02 eV. We use an anisotropic 



cutoff energy for the local field effects . Since the sur- 
face plasmon depends sensitively on the density profile at 
the surface where the density decays exponentially into 
the vacuum, a cutoff energy of 500 eV is applied in the 
z-direction. Compared to the RPA results in panel (a), 
the inclusions of the LDA exchange-correlation kernel in 
panel (b) shifts the peaks down by 0.1 - 0.2 eV. 

The energies of these surface plasmons for both the 
T — M and T — K directions are shown in Fig. [5] (c). 
The obtained dispersion relations agree well with previ- 
ous calculations^. The well known negative dispersion 
at small q observed for simple metal surfaces are also 
well reproduced in this work. Compared to experimental 
data, the TDLDA energies of the surface plasmons agree 
within 0.1 eV for small q, while the discrepancy increases 
to around 0.2 eV for larger q, which can attributed to the 
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FIG. 7: (Color online) Loss function of free standing graphene 
(a) and graphene on SiC substrate (b) as a function of q. The 
loss functions, from bottom to top (solid lines), correspond to 
increasing q at an interval of 0.046/ A. The dashed line corre- 
sponds to the loss function of the substrate at g = 0.092/ A. 
(c) Dispersion relations for the vr plasmons of free standing 
graphene (red filled circles) and graphene on SiC (black filled 
squares). They are compared with earlier ab-initio calcu- 
lation on free standing graphene (blue hollow circles) and 
experiments on single wall carbon nanotubes (green hollow 
squares)^ as well as experiments on graphene / SiC(OOOl) 
(purple hollow diamonds)^^. Lines are added to guide the 
eye. 



FIG. 6: (Color online) Atomic structure of graphene adsorbed 
on SiC(OOOl) (a-fb) and Al(lll) (d-|-e). The lateral unit cells 
are indicated by red lines in the top panels. The LDA band 
structures of the surfaces are shown in the lower panels. Also 
shown is the band structure of free standing graphene (red 
dots). The Fermi level is set to zero. 

fact that the ALDA kernel is q-independenl4i. 

D. Plasmons in adsorbed graphene 

In this section we investigate the influence of a sub- 
strate on the plasmon excitations in graphene. For a 
more detailed discussion of these results we refer the 
reader to Ref. Hsj. As representatives for semiconduct- 
ing and metallic substrates we consider SiC (0001) and 
Al(lll). Both of these systems bind the graphene rela- 
tively weakly so that hybridization effects are relatively 
unimportant. Thus the largest effect of the substrate is 
expected to arise from the long range Coulomb interac- 
tion between electrons in the two subsystems. 

The atomic structure and band structure of graphene 
on both substrates are shown in Fig. [51 For 
graphene/SiC(0001), the unit cell, indicated by red sohd 
lines in panel (a), contains 2x2 graphene and \/3 x \/3 
SiC— As can be seen in panel (b), two carbon layers 
are adsorbed on four bi-layers of SiC and the dangling 



bonds at the backside of the slab are saturated by hydro- 
gen. The first carbon layer adsorbs covalently on the SiC 
surface and is here considered as a part of the substrate. 
The upper carbon layer binds weakly to the substrate, in 
agreement with experiments^, with an LDA binding en- 
ergy per C atom of 0.039 eV, and adsorption distance of 
3.56 A. As shown in panel (c), linear conical bands appear 
within the bandgap of the substrate, resembling that of 
free-standing graphene (red dotted line). The Fermi level 
is shifted up by 0.05 eV, introducing slight electron dop- 
ing into graphene. For the graphene/ Al( 111) structure 
we use a 1 X 1 unit cell with four layers of Al as sub- 
strate. Again, graphene binds weakly to the Al surface 
with an LDA interplane distance of 3.36A and binding 
energy per C atom of 0.049 eV, in good agreement with 
recent van der Waals DFT calculations^^. As shown in 
panel (f), the 'Dirac cone' of graphene is shifted 0.5 eV 
below the Fermi level. The computational details for cal- 
culation of the loss functions can be found in Ref. |45| . 

Fig. [7] shows the calculated loss function of free stand- 
ing graphene (a) and graphene on SiC (b). The free- 
standing graphene exhibits a collective mode at around 
5 eV, which results from the electronic transitions of the 
TT — > TT* bands and is referred to as the graphene tt plas- 
mon. The dispersion of the tt plasmon is shown in Fig. 
[Tfc). In contrast to its three dimensional counterpart, 
graphite, which shows a parabolic dispersion of the tt 
plasmons^"^, graphene has a linear plasmon dispersion. 
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FIG. 8: (Color online) Surface loss functions for graphene 
on Al(lll) as a function of the adsorption distance d for 
|q| = 0.046/A. The surface loss function of free standing 
graphene is shown as black dots. Inset: sketch of graphene 
on Al substrate. 



The origin of the linear dispersion has been attributed to 
the role of local field effects^. 

Fig. [7] (b) shows the loss function of graphene ad- 
sorbed on the SiC(OOOl) surface. Compared to the re- 
sults of the free standing graphene, the strength of the 
TT plasmons are strongly damped, in particular for small 
q values. As q increases, the strength of the tt plasmons 
gradually recovers to that of a free standing graphene, 
indicating that the substrate effect becomes weaker for 
larger q. As shown in Fig. [7|c), the substrate has little 
effect on the energies of the tt plasmons. In fact the plas- 
mon dispersion for both free standing and substrate sup- 
ported graphene agree well with previous calculations^^ 
as well as experiments on graphene/SiC^ and carbon 
nanotubes^. We have found that the response func- 
tion, and thus the EELS spectrum, of the combined 
graphene/substrate system can be obtained accurately 
from the response functions of isolated graphene and sub- 
strate assuming only Coulomb interaction between the 
two, i.e. neglecting effects related to hybridization and 
charge transfer—. This demonstrates that the strong 
damping of plasmons results from the non-local screen- 
ing of the graphene plasmon excitation by the substrate 
electrons. 

Fig. [5] shows the surface loss function of graphene on 
Al(lll) for various adsorption distances. In contrast to 
the semiconducting SiC substrate, the tt plasmon at 5 
eV is completely quenched on the metallic Al substrate 
at the equilibrium distance d — 3.36 A(full black line). 
As the graphene is pulled away from the surface, the 



TT plasmon reappears at an energy lower than that of 
the free standing graphene. This downshift is due to 
the coupling to the surface plasmons of the aluminum 
substrate at 9.0 eV. The graphene tt is fully recovered 
at a distance of around 20 A illustrating the long range 
nature of the interaction. 



V. CONCLUSIONS 

We have implemented the linear density response 
function in the adiabatic local density approximation 
(ALDA) within the real space projector augmented wave 
method GPAW, and used it to calculate optical and di- 
electric properties of a range of solids, surfaces and in- 
terfaces. The Kohn-Sham wave functions, from which 
the response function is built, can be obtained either on 
a real space grid or in terms of localized atomic orbital 
basis functions. The latter option reduces the compu- 
tational requirements for calculating and storing the of- 
ten very large number of wave functions required for the 
construction of the response function without sacrificing 
accuracy. The dielectric constants of a number of bulk 
semiconductors as well as the optical absorption spec- 
trum of silicon at the ALDA level was shown to be in 
good agreement with previous calculations. For the sur- 
face plasmons of the Mg(OOOl) surface we find, in agree- 
ment with previous studies, that the ALDA kernel low- 
ers the plasmon energies by around 0.3 eV realtive to 
the RPA values and thereby reduces the deviation from 
experiments from around 4% to 1-2%. Very good agree- 
ment with experiments was also found for the plasmon 
energies of graphene which were shown to exhibit a lin- 
ear dispersion with a value of 4.9 eV in the long wave 
length limit. The deposition of graphene on a SiC sub- 
strate is shown to have little effects on the plasmon en- 
ergies but leads to significant damping of the plasmon 
resonances. In contrast deposition on an Al surface com- 
pletely quenches the graphene plasmons due to strong 
non-local electronic screening effects. 
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